#### Produce regression tables for main and alternative spec.

### Prepare environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- F
do.top100 <- F

source("load.R", encoding="iso-8859-1")

library(lfe)
library(splines)

if (do.top100) {
    library(dplyr)
    top100 <- (df %>% group_by(Munic�pio) %>% summarize(total.produce=sum(total.produce)) %>% arrange(-total.produce))$Munic�pio[1:100]
    subdf <- subset(df, Munic�pio %in% top100)
} else {
    subdf <- df
}

## Evaluate temporal and spatial correlation distance
disteval <- cbind(subdf[, c("Munic�pio", "year", "latitude", "longitude")], resid=NA)
disteval$resid[-mod6$na.action] <- mod6$resid
deres.lag <- data.frame()
for (muni in unique(subdf$Munic�pio)) {
    rows <- which(disteval$Munic�pio == muni & !is.na(disteval$resid))
    if (length(rows) > 1) {
        acf_results <- acf(disteval$resid[rows], plot = FALSE)
        sig_acf_lags <- which(acf_results$acf > 1.96/sqrt(length(rows)))
        deres.lag <- rbind(deres.lag, data.frame(Munic�pio=muni, lags=length(sig_acf_lags)))
    }
}

library(gstat)
library(rgdal)
deres.dis <- data.frame()
for (year in unique(subdf$year)) {
    rows <- which(disteval$year == year & !is.na(disteval$resid) & !is.na(disteval$longitude))
    subde <- disteval[rows,]
    if (length(rows) > 1) {
        coordinates(subde) <- ~ longitude+latitude
        proj4string(subde) <- CRS("+proj=longlat +datum=WGS84")
        subde2 <- spTransform(subde, CRS("+proj=utm +zone=22 +datum=WGS84"))
        varfit <- variogram(resid ~ 1, subde2)
        varfit2 <- fit.variogram(varfit, vgm("Sph"))
        deres.dis <- rbind(deres.dis, data.frame(year, range=varfit2$range[2] / 1000))
    }
}

## Run the regressions

if (do.origspec) {
    mod1 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod2 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod3 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod4 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod5 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod6 <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod1b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio | 0 | Munic�pio, data=subdf)
    mod2b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
    mod3b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio | 0 | Munic�pio + year, data=subdf)
    mod4b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
    mod5b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio | 0 | Munic�pio + year, data=subdf)
    mod6b <- felm(logyield ~ tmin + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
} else {
    mod1 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod2 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod3 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod4 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod5 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio | 0 | Munic�pio, data=subdf, keepCX=T)
    mod6 <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio + year | 0 | Munic�pio, data=subdf, keepCX=T)
    mod1b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio | 0 | Munic�pio, data=subdf)
    mod2b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
    mod3b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio | 0 | Munic�pio + year, data=subdf)
    mod4b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:year | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
    mod5b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio | 0 | Munic�pio + year, data=subdf)
    mod6b <- felm(logyield ~ tmin + tmax + gdd1000 + edd1000 + prcp + prcp2 + state:ns(year, df=3) | Munic�pio + year | 0 | Munic�pio + year, data=subdf)
}

library(stargazer)

## Write these to yieldmod-3se[-top]-SUFFIX
## Stars according to first model (b - M & Y Cl.)

stargazer(list(mod1b, mod2b, mod3b, mod4b, mod5b, mod6b), add.lines=list(c('Municipality FE', rep('Yes', 6)),
                                                       c('Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes'),
                                                       c('State Trends', 'None', 'None', 'Linear', 'Linear', 'Cubic', 'Cubic')),
          dep.var.labels="Log Yields", covariate.labels=weatherlabels,
          omit=':', df=F, float=F)

## Place SEs above
stargazer(list(mod1, mod2, mod3, mod4, mod5, mod6), add.lines=list(c('Municipality FE', rep('Yes', 6)),
                                                       c('Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes'),
                                                       c('State Trends', 'None', 'None', 'Linear', 'Linear', 'Cubic', 'Cubic')),
          dep.var.labels="Log Yields", covariate.labels=weatherlabels,
          omit=':', df=F, float=F)

source("conley.R")

cse1 <- ConleySEs(mod1, subdf[-mod1$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
cse2 <- ConleySEs(mod2, subdf[-mod2$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
cse3 <- ConleySEs(mod3, subdf[-mod3$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
cse4 <- ConleySEs(mod4, subdf[-mod4$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
cse5 <- ConleySEs(mod5, subdf[-mod5$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
cse6 <- ConleySEs(mod6, subdf[-mod6$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)

stargazer(list(mod1b, mod2b, mod3b, mod4b, mod5b, mod6b), add.lines=list(c('Municipality FE', rep('Yes', 6)),
                                                       c('Year FE', 'No', 'Yes', 'No', 'Yes', 'No', 'Yes'),
                                                       c('State Trends', 'None', 'None', 'Linear', 'Linear', 'Cubic', 'Cubic')),
          dep.var.labels="Log Yields", covariate.labels=weatherlabels,
          omit=':', df=F, se=list(sqrt(diag(cse1$Spatial_HAC)), sqrt(diag(cse2$Spatial_HAC)), sqrt(diag(cse3$Spatial_HAC)), sqrt(diag(cse4$Spatial_HAC)),
                  sqrt(diag(cse5$Spatial_HAC)), sqrt(diag(cse6$Spatial_HAC))))
